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ABSTRACT 

Reverberation mapping probes the structure of the broad emission-line region (BLR) 
in active galactic nuclei (AGN). The kinematics of the BLR gas can be used to measure 
the mass of the central supermassive black hole. The main uncertainty affecting black 
hole mass determinations is the structure of the BLR. We present a new method 
for reverberation mapping based on regularized linear inversion (RLI) that includes 
modelling of the AGN continuum light curves. This enables fast calculation of velocity- 
resolved response maps to constrain BLR structure. RLI allows for negative response, 
such as when some areas of the BLR respond in inverse proportion to a change in 
ionizing continuum luminosity. We present time delays, integrated response functions, 
and velocity-delay maps for the H /? broad emission line in five nearby AGN, as well as 
for Ha and Hy in Arp 151, using data from the Lick AGN Monitoring Project 2008. 
We find indications of prompt response in three of the objects (Arp 151, NGC 5548 
and SBS 1116+583A) with additional prompt response in the red wing of H /?. In SBS 
1116+583A we find evidence for a multimodal broad prompt response followed by a 
second narrow response at 10 days. We find no clear indications of negative response. 
The results are complementary to, and consistent with, other methods such as cross 
correlation, maximum entropy and dynamical modelling. Regularized linear inversion 
with continuum light curve modelling provides a fast, complementary method for 
velocity-resolved reverberation mapping and is suitable for use on large datasets. 
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1 INTRODUCTION 

Active galactic nuclei (AGN) are believed to be powered by 
the release of gravitational potential energy when matter 
falls onto supermassive black holes in the centres of galax¬ 
ies. Some AGN have broad emission lines that are thought to 
be Doppler broadened emission from gas orbiting the central 
black hole. The broad lines vary in response to the contin¬ 
uum emission, suggesting that they are powered by ionizing 
radiation originating in the immediate vicinity of the black 
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hole ([Bahcall &; Kozlovskvlll969l : lDavidsonlll972h . The time 
delay between the variations in the continuum and emis¬ 
sion lines can be used to measure the structure and char¬ 
acteristic size of the broad emission li ne region (BLR) by 
the method of reverberat i on mapping ([Bahcall et alJ 1 19721 : 
iBlandford fe McKeelll983 : |Petersonlll993l ) . 


The radius of the BLR combined with the width of 
the emission line provi des a measurement of the mass of 
the c entral black hole ([Peterson et alJ Il998bl : iKasoi et al.l 
l200(tl . Reverberation masses have been found to corre¬ 
late with properties of the host galaxy, su ch as stel¬ 
lar velocit y dispersion (e.g. WoojE+jd. 2010) and bulge 
mass (e.g. iBentz et al.l l20Q9al : iKormendv fe Hoi 1120131 ) , as 
has been found for black holes in inactive galaxies 
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Rchstone 19951: iMagorrian et al. 1998: 
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2000 : iGebhardt et al.l l200olj and indi- 


eating a connection between supermassive black hole growth 
and galaxy evolution. In addition, BLR radii have been 
found to correlate with the UV/optical luminosity of th e 
AGN (e.g. lKasni et~aD200(ll20nalBent,z et, al.l20n9bl . Efn^ . 
suggesting that AGN can be used as independent stan¬ 
dard candles for cosmology ([Watson et al.ll201ll : iKing et al.l 
l2014h . The radius-luminosity relation also enables black hole 
masses to be estimated out to large redshift (z > 6) by mea¬ 
suring the AGN luminosity and broad emission line width 
in a single spec trum, yielding sing l e-epoch black hole mass 
estimates (e.g. IVe stergaard 2002; McL ure &; Jarvisl 120021 : 
IVestergaard fe Peterson|[2006: IShenll2013lb 

Because of limited knowledge of the structure of the 
BLR, it is necessary to introduce a scaling factor in tra¬ 
ditional reverberation mapping black hole mass measure¬ 
ments, such that the black hole mass is determined by 


Mbh 


f 


AV^blr 

G 


( 1 ) 


where AV is the velocity width of the varying part of the 
emission line, G is the gravitational constant, and / is the 
order-unity scaling factor that encapsulates th e unknown de¬ 
tails r egarding the BLR gas and kinematics (W andel et al.l 
Il999l : lOnken et al.l l2004h . The /-factor is generally cali¬ 
brated using the Mbh — relation between the black hole 
mass and the stellar velocity dispersion of the bulge in the 
host galaxy. This relies on determining Mbh independently 
in quiescent galaxies using stellar and gas dynamics (see 
iKormendv &; Holl2013l for a review), and then determining 
the average scaling factor (/) that brings the ensemble of 
active galaxies into agreement with the quiescent Mbh — <r* 
relation. This process only corrects for the average offset, 
but it results in black hole masses that are uncertain for 
any particular AGN by a factor of 2 — 3. In addition, (/) 
is determined using H /3, but often applied to other emis¬ 
sion lines such as Mgn A2799 and Civ A1549, although the 
structur e of the BLR for th e se emission lines is very un¬ 
certain (|Metzroth et al]l2Q0(i : IWooll2008l : iKasoi et "all [20071 : 
iTrevese et al JI 201 J. 

The main uncertainty in measuring black hole masses 
using reverberation mapping thus comes from our limited 
knowledge of the kinematics and geometry of the BLR. A 
precise measurement of the kinematic and geometric struc¬ 
ture of the BLR would enable the determination of / for 
any individual AGN. This could help reduce scatter in black 
hole mass determinatio ns for all single epoc h mass estimates 
and their applications (| Kelly Sz Shenll2013h . A better under¬ 
standing of changes in BLR structure with luminosity and 
redshift will also reduce systematic errors possibly affecting 
current scaling relations. To make a full map of the BLR we 
need to determine not just its characteristic radius, but also 
the full velocity-resolved transfer function that describes the 
relation between the continuum emission and the emission 
line response. 

To first order, the problem of reverberation mapping 
can be formulated as a deconvolution problem in which the 
flux in the emission-line light curve, at a given wavelength A, 
Fi (t, A) is given by a convolution of the AGN continuum light 
curve, over some wavelength range, F c (t) with a transfer 
function T(t, A) that encodes the physics and geometry of 


the BLR, 


/ oo 

*(t, X)F c (t — t) dr. 

-OO 


( 2 ) 


The transfer function, as a function of time delay and 
wavelength T(y , A), is called the velocity-delay map 
( We lsh V Hornel Il99lh . In this approximation, the trans¬ 
fer function is assumed to be linear. Even though the de¬ 
tailed physics is likely to be more complicated, the linear ap¬ 
proximation is currently sufficient given that observational 
datasets have only recently become good enough to probe 
the full velocity-resolved transfer function. Thus, the trans¬ 
fer function introduced here represents an observed projec¬ 
tion of the underlying structure of the AGN. A sound in¬ 
ference of the transfer function will need to account for any 
residual mismatches between the assumed model and the 
data, such as non-linearities. 

Because transfer functions represent projections of 
the underlying physical structure, physical and geometri¬ 
cal models are required to interpret them. Several groups 
have gone through the exercise of predicting transfer func¬ 
tions based o n underlying physical models for the BLR 
structure (e.g. Chiang fe Murray Jl996l : iBottorlf et al.|[l997l : 
iGoad et alll2012l : IPerez et"^L 19921 ). These studies provide 
a valuable catalogue of transfer functions that can be con¬ 
sulted when interpreting results from reverberation map¬ 
ping. 

Much effort has gone into developing efficient methods 
for estimating the BLR size and transfer function T (r, A). 
Early attempts relied on estima ting the time delay by-eye 
(IChereoashchuk &; Lvutvil Il973h , but many sop h isticated 
metho ds have been developed since. iBlandford &; McKeel 
(Il982l j were the first to calculate transfer functions directly 
using the convolution theorem of Fourier transforms and 
provided a comprehensive catalogue of transfer functions for 
a number of idealised BLR structures. Unfortunately, the re¬ 
quirement for very high quality data, as well as difficulties 
in dealing with measurements errors, mean that the method 
has seen little application since its publication. This may 
change with future high cadence, high signal-to-noise rever¬ 
beration mapping campaigns. 

Another type of inversion procedure is the maximum 
entropy method that finds the solution for the transfer func¬ 
tion that has the hig hest entropy, while stil l providing a 
good fit to the data (ISkilling &; Brvar] 1 19841 ; IKrolik et al.l 
Il99ll : iHorne et al.l Il991 : Ih ornel Il994| j. Maximum entropy 
has been successful in estimating transfer f unctions and 
veloc i ty-delay maps in a number of AGN (e.g. IKrolik et al 
Ulrich &; Ho rne 19961: iBentz et al.l 1 201 Obi : I Grier et aT 


199ll: Jl 

2013af r 


2013al j. The downsides of maximum entropy are that it is 


computationally expensive, it relies on a number of assump¬ 
tions about the shape of the transfer function, and it is dif¬ 
ficult to carry out rigorous error analysis and model com¬ 
parisons. 

Another class of method is dynamical modelling in 
which a full physical model of the BLR is constructed, and 
its parameters are inferred within the framework of Bayesian 
statistics. The statistical framework allows for rigorous er - 
ror analysis as well model selection (IPancoast et al.l [201 lh . 
Furthermore, dynamical modelling circumvents the need for 
interpreting transfer functions by providing direct estimates 
of physical model parameters such as inclination and the 
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black hole mass, which in turn allows for the /-factor to be 
calculated directly. The main challenges of dynamical mod¬ 
elling methods are that they require long computation times 
and the assumption that the model is flexible enough to pro¬ 
vide a good description of the BLR. 

Dynamical modelling and maximum entropy both pro¬ 
vide useful and complementary constraints on the BLR 
structure. They also rely on a number of assumptions about 
the allowed shape of the transfer function and are fairly 
computationally expensive. It is therefore worthwhile to con¬ 
sider alternative methods for analysing reverberation map¬ 
ping data that are less computationally expensive and allow 
for greater flexibility in estimating velocity-resolved transfer 
functions. 

Here we develop a method for reverb eration ma pping 
based on regularized linear inversion (RLI: Ivio et al.l Il994l : 
iKrolik &; Donelll995l ) , which we extend by including statisti¬ 
cal modelling of the AGN continuum light curve light curves. 
The method is complementary to other reverberation map¬ 
ping techniques, and has the advantage that it provides an 
analytical expression for the transfer function with very few 
assumptions about its shape. RLI is a flexible, free-form 
method, allowing the data to suggest the form of the inferred 
transfer function. This makes RLI an ideal tool for exploring 
BLR physics beyond the framework of current BLR models. 
At the same time RLI is one of the fastest reverberation 
methods that provides a direct estimate of the transfer func¬ 
tion. We extend the RLI method to include error in the light 
curves by a combination of Gaussian process modelling and 
bootstrap resampling, thereby providing a robust estimate 
of the transfer function and its uncertainties. 


As a first application, we apply our method to photo¬ 
metric and spectroscopic light curves of five nearby AGN 
measured by the Lick A GN Monitoring P roject 2008 col¬ 
laboration (LAMP 2008: IBentz et al1l2008l ). The main pur¬ 
pose of this project was to measure masses of supermas- 
sive black holes in 13 nearby (z < 0.05) Seyfert 1 galaxies 
(jBentz et alJl2009ch . Besides improving black hole mass es¬ 
timates, LAMP 2008 has produced a medley of scientific re¬ 
sults including: AGN variability characteristics (| Walsh et al.l 
l2009l b an update of t he Mbh — cr* rel ation with reverber¬ 
ation mapped AGN (IWoo et alJl20ld), de tailed reverber¬ 
ation mapping studies (jBentz et* 'al.l l2010alfb ). prob ing the 
Rblr — L relation in the X-rays ( Greene et alJl2010h. recal¬ 
ibrating single-epoch black hole masses (IPark et al. 20121) 
and dy namical modelling of the H j3 BLR ([Pancoast et al.l 
l2014bh . 


We analyse five objects from LAMP 2008: Arp 151 (Mrk 
40), Mrk 1310, NGC 5548, NGC 6814 and SBS 1116+583A, 
providing integrated response functions, time delays and 
velocity-delay response maps for the H /3 emission line in 
each object, as well as Ha and Hy in Arp 151. We be¬ 
gin in Section [2] by describing how we obtained light curves 
from the LAMP 2008 dataset. In Section [3] we outline the 
RLI method and show results from simulations. Section 0] 
presents the main results of our analysis. The results are 
discussed in Section [5] and finally we provide a short sum¬ 
mary and conclusions in Section [6] 


All time delay-axes in the figures are in the observed 
frame. All time delays quoted in the text and in Table 1 are 
in the rest frame of the AGN. 


2 DATA 

We use data from LAMP 200^3, a dedicated spectroscopic 
reverberation mapping campaign that ran for 64 nights at 
the Lick Observatory 3-m Shane telescope. The spectro¬ 
scopic data were supplemented by photometric monitoring 
using smaller telescopes, including the 30-inch Katzman Au¬ 
tomatic Imaging Telescope (KAIT), the 2-m Multicolor Ac¬ 
tive Galactic Nuclei Monitoring telescope, the Palomar 60- 
inch telescope, and the 32-inch Tenagra II telescope. A de¬ 
tailed description of the LAMP 2 008 observing camp aign 
and initial results a re published in Walsh et al.l (|2009h and 
iBentz et al.l (|2009ch . 

2.1 Continuum light curves 

For our analysis, we use Johnson B and V broad band 
continuum light curves from LAMP 2008. The bands were 
chosen to improve dynamical modelling results and resolve 
variability features. The fl uxes are measured using standard 
aperture photometry (see I Walsh et al.l 120091 for a complete 
discussion). The light c urves are the same as t hose used in 
dynamical modelling bv IPancoast et all (l2014bl ). The B and 
V b and light curves ar e very similar for all objects analysed 
(see IWalsh et al.l 12009 ), and the choice of continuum light 
should not significantly affect our results. 


2.2 Emission line spectra and light curves 

Emission-line light curves are generated from flux-calibrated 
spectra from the LAMP 2008 campaign. The H [3 emission 
line is isolated in all objects using spectral decomposition, 
by modelling all line and continuum components individu¬ 
ally, and subtr acting away e verything but the emission line 
of interest (|Park et alJl2012h . We decide to keep the narrow 
component of H [3 to avoid introducing additional error at 
the centre of the line. This should not affect our results, as 
the narrow-line component i s constant o n the time-scale of 
the observing campaign (see IPark et alJl2Q12h . and because 
we only consider variations around the mean flux in the line. 
In Arp 151 we also analyse the Ha and Hy lines, where Hy 
is isolated usi ng spectral decomp osition in a way similar to 
that for H /3 dBarth et al.l 1 20 111 ) . Because no spectral de¬ 
composition is avai lable for Ha we follow the procedure of 
IBentz et al.l (2010a ). isolating Ha by subtracting a straight 
line fit to two wavelength windows on either side of the line. 

The resulting emission line spectra have the same spec¬ 
tral resolution as the origina l LAMP 2008 data. T he spectral 
dispersions are provided in IBentz et al.l ( 2009c ) and range 
from 11.6 — 14.7 A (FWHM), corresponding to 5.9 — 7.5 pix¬ 
els per resolution element. 

A few of the spectra in the LAMP 2008 dataset have 
been identified as unreliable due to l ow spectral qual ity or 
suspicious features above the noise ([Park et al.ll2012l ). We 
tested the effect of removing these spectra in the analysis 
and found that it made little difference to the overall results. 
Because of this, and to be able to compar e directly with 
recent results from dynamical modelling ([Pancoast et al .1 
l2014bh . we decided to include the unreliable spectra in the 
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Data available at http: //www. physics. uci. edu/~barth/lamp. html 
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RLI analysis presented here. The emission line light curves 
for H [3 are thus ide ntical to those used for the dynamical 
modelling analysis in lPancoast et all dl2014b j. 


3 REGULARIZED LINEAR INVERSION 

Regularized linear inversion seeks to solve the transfer equa¬ 
tion (Equation [2]) for the transfer function T(r, A) analyti¬ 
cally and without any assumptions regarding its functional 
form (only that it is a bounded linear operator). This ap¬ 
proach is potentially very powerful in that it relies solely 
on the data when deriving the transfer function. This in 
turn allows for very little freedom for the method to select 
a proper solution. In the presence of noise, or if the system 
deviates slightly from the linear assumption, this lack of free¬ 
dom means that a solution cannot be found at all. Therefore, 
instead of looking for a unique solution that fits the data, 
the method determines the solution that minimizes the x 2 - 
A simple x 2 minimization will generally over-fit the data and 
make the inversion unstable. This is overcome by a regular¬ 
ization in which the first-order derivative of the solution is 
to be minimized along with the y 2 , such that the solution is 
smoo thed at the level of the noise (|Phillipsl [l962; Tikhonov) 
Il963h . The result is comparable to maximizing the entropy 
under a y 2 constraint as in the case of maximum entropy 
methods. RLI has several advantages over other methods, 
specifically it 1) makes no assumption about the shape or 
positivity of the transfer function, 2) can be solved analyti¬ 
cally, and 3) has very few free parameters (the regularization 
scale as well as the transfer function window and resolution). 
Note that while RLI can in principle fit any shape of trans¬ 
fer function, in practice the result will be limited by the 
sampling and noise in the data, which means that a mini¬ 
mum scale exists below which the response function will be 
unresolved. 

By measuring the continuum (F c ) and emission line 
(ifi) light curves we can in principle solve the transfer equa¬ 
tion m for T(t, A). In reality, data is always discrete, so we 
rewrite the transfer equation as a linear matrix equation of 
the form 


about the structure of the BLR relies on a number of as¬ 
sumptions. First, we assume that the variations in the AGN 
continuum bands are correlated with the AGN ionizing con¬ 
tinuum (this is not necessarily true, see e.g. ICollier et all 
Il998l : iMcHardv et al.ll2014l : lEdelson et al.ll2015l who find ev¬ 
idence for a time delay between the UV and optical contin¬ 
uum in NGC 7469 and NGC 5548). Second, we assume that 
the continuum emission originates from a region negligible 
in size compared to the BLR. Third, we assume that the 
continuum ionizing radiation is emitted isotropically. Last, 
we assume that the BLR structure is constant and the re¬ 
sponse linear for the duration of the campaign, such that 
a single linear transfer function can be calculated from the 
full light curves. 


3.1 Solving for the response function 


Rather than solving for the transfer function, which includes 
constant emissivity components of the BLR, we consider 
only variations about the mean of the light curves. By doing 
this, we are measuring the response of the line emission to a 
change in the continuum emission that is ionizing the BLR 
gas. Hence the qu antity we are solv i ng for, T(r, A), is the 
response function (|Krolik et al.lll99ll: iGoad et a 1ll993h . 

Following the notation of IKrolik &; PoneT i 1995b and 
considering only variation about the mean of the light 
curves, we write x 2 as 


N 1 


SFi(U) ~ J21 Fc(U - Tj ) - (F c )]'b(tj) 

3 = 1 


This expression can be recast to matrix notation, 


(4) 


X 2 = (L - C*) 2 , (5) 


where the light curves enter as 


C ij = [F c (t-T j )-{F c )]/ai(U) (6) 

L = SFi(U)/m(ti), (7) 


and the variation about the mean in the emission-line light 
curve is defined as 


Laa = ^aaC, (3) 

where Laa is the emission line light curve integrated over the 
wavelength range (spectral bin) A A, C is a matrix of contin¬ 
uum light curves (see below), and ^aa is the transfer func¬ 
tion corresponding to the given wavelength range. While we 
allow the transfer function to depend on wavelength, sym¬ 
metries in the BLR as well as observational projections will 
tend to correlate transfer functions in the time and wave¬ 
length domains. The set of transfer equations over a range 
of frequencies across an emission-lin e we call the velocity - 
delay map for the given emission line (IWelsh fc Hornel[l99lh . 

For perfect noise-free data solving the discrete trans¬ 
fer equation © would be a simple matter of inverting C 
to obtain the transfer function T r . Because of noise we can¬ 
not hope to find an exact solution to the linear inversion 
problem. Instead we seek a solution that minimizes the x 2 
together with a smoothing condition that ensures that we 
are not fitting the noise. 

Interpreting the transfer function to make statements 


8F l (t i ) = F l (ti)-(Fi). (8) 

We choose a time delay resolution of 1 day, as this corre¬ 
sponds to the highest resolution of the data. By running the 
analysis with different resolutions we have confirmed that 
the choice of resolution, within reasonable value s , doe s not 
change our results. Contrary to IKrolik &; Don el |l995), we 
calculate a constant mean of the continuum light curve data 
points such that (F c ) is not a function of the time delay 
r. This is done in order to allow meaningful comparisons 
between different continuum realizations (see Section 15311 . 
Minimizing \ 2 in Equation © leads to the linear equation 

C T C^ = C t L. (9) 

Although this expression looks simple, it turns out to be 
ill-conditioned. To remedy this, we put an extra constraint 
on the problem, namely that the solution should be smooth 
at the scale of the noise (this effectively avoids fitting the 
noise). To guarantee smoothness of the solution, we intro¬ 
duce a differencing operator H acting on \l>, and require that 

































Constraints on the BLR from regularized inversion 5 


the first-order difference (the discrete version of the first- 
order differential) be minimized together with the ordinary 
X 2 . To control the weights between the x 2 and the first-order 
difference, a scaling parameter k is introduced that sets the 
scale of the regularization. Thus, the expression to minimize 
becomes 


(C T C + kH t H)$ = C t L. (10) 


This expression is more stable under inversion, while sac¬ 
rificing detail by emphasizing smoothness of the solution. 
The question is then how to choose a suitable scale Hi for the 
regularization. 

The best choice of regularization scale depends on the 
signal-to-noise in the data, as well as the level of uncorre¬ 
lated systematic uncertainties that result in deviations from 
the assumption of linear response. As a starting point for 
selecting a regu l arizat ion scale k, we follow the recommen¬ 
dation of lPres 3 CuH) in which k, — Ko is chosen to provide 
equal weights to the two left-hand-side terms in equation 

m, 


_ Tr(C T C) 
" Tr(H T H)' 


( 11 ) 


iKrolik Sz Done! ([l995) suggest using the larges t value of k 
which gives an acceptable x 2 > while IVio et all (|l994h sug¬ 
gest using a value of Hi which provides the “best” trade-off 
between resolution and noise. By running tests on simulated 
data, we find that good results are generally achieved for 
k = k>o. For this reason, and to reduce the number of free 
parameters in the method, we fix k = Ko for all our anal¬ 
ysis (see Section 13.61 for a test of the effect of changing the 
regularization scale). 

We calculate all response functions in the time delay 
interval 0 — 30 days. The lower bound is chosen to impose 
causality in the solution. The upper bound is chosen to be 
well beyond the time-scale where we expect a significant re¬ 
sponse based on pr evious measureme nts and expectations 
for Seyfert galaxies ([Bentz et al.ll2013h . We tested the effect 
of changing the time delay windows and found that it had 
only minor effects on the results, as longs as the window is 
long enough to include the main response. When presenting 
the results, we show only the first part of the response func¬ 
tion from 0 — 15 days, as we found the significant response 
power to be located at these scales. 


3.2 Continuum light curve errors 

Continuum flux errors are not included explicitly in the 
RL1 formalism presented above. Instead, we include con¬ 
tinuum errors and interpolation by modelling the contin- 
uum light curves using Gaussian processes ([Pancoast et alJ 
B)- Studies of AGN variability, using sampling intervals of 
days, have suggested that AGN continuum variability is well 
modelled by a da mped random wal k or Ornstein-Uhlenbeck 
(O-U) process (iKellv et al.l I2QQ9I: iKozlowski et al] l20ld : 
iMacLeod et~aD l20ld : IZu et al.l l2013h . Formally this is a 
CAR(l) or continuous-time first-order autoregressive pro¬ 
cess, which is a stationary Gaussian process with a power 
spectral density (PSD) slope of —2. The CAR(l) process is 
often used in reverberation mapping methods for modelling 
the continuum light curve and enable efficient interpola¬ 


tion and error estimation (e.g. IZu et al.ll201ll : IPancoast et al.l 

B)- 

Recent hig h cadence observations of AGN using the Ke¬ 
pler spacecraft (|Mushotz kv et al .ll201ll:ICarini &; Rvldl2012l : 
lEdelson et, a. 1.1 [20141 : iKa.sliwal et, al .11201 oil . have found steep 

PSD slopes (< —2) inconsistent with a CAR(l) process. The 
effect of having a steeper PSD slope is to suppress small scale 
variability in the light curves. Using a CAR(l) process for 
modelling the continuum light curves may therefore result 
in over-fitting spurious features, due to noise at short time 
scales. This introduces artificial structure at the level of the 
noise, that can lead to underestimated errors being prop¬ 
agated form the continuum model to other reverberation 
mapping parameters. 

The magnitude of the error introduced by assuming the 
wrong continuum model depends on the sampling rate and 
errors in the data, as well as how it is implemented in the re¬ 
verberation mapping method. Regularized linear inversion, 
as implemented here, mitigates the effects noise in the light 
curves by smoothing the solution at the shortest time scales. 
This results in a solution that is dominated by the longer 
multi-day variability features. For this reason, the response 
functions derived should not be affected by the details of the 
assumed continuum model, as long as multi-day time scale 
features are accurately reproduced. Even so, the shallower 
slope of the CAR(l) process, compared to Kepler AGN light 
curves, may still lead to underestimated errors. 

To determine how the choice of continuum model PSD 
affects our results, we run RLI using a number of dif¬ 
ferent continuum models with varying degrees of small 
scale structure, corresponding to varying PSD slopes. We 
do this by changing the power a in the covariance func¬ 
tion for the Gaussian process. The covariance is given by 
([Pancoast et al.ll201 lh 


C(ti,t2) = cr 2 exp 


|t 2 - ti\ 


( 12 ) 


where a is the long-term standard deviation, r is the typical 
time scale of variations and the power a take on values in the 
interval [1, 2]. A power of a = 1 corresponds to the CAR(l) 
model, whereas a > 1 produces smoother continuum mod¬ 
els, with PSD slopes < — 2. Using data from LAMP 2008 for 
Arp 151 (see results in section [4. 2. ID . we test four different 
continuum models with a = 1.0, 1.8, 1.9 and 2.0 respec¬ 
tively (see Figure [I]). We find that the emission line fits and 
response functions are not strongly affected by the assumed 
continuum model. For this reason, and to be consistent with 
other reverberation mapping methods such as javelin and 
dynamical modelling, we choose to use a CAR(l) process as 
a prior for the Gaussian process continuum modelling. Fu¬ 
ture analyses of high-cadence reverberation mapping data 
should likewise be careful in considering the effects of the 
assumed continuum model. 

The finite sampling and duration of the light curves 
hampers our ability to model structure on scales close to 
or below the sampling time scale, and on scales longer than 
the duration of the light curve. At short time scales the light 
curves will likely be dominated by observational errors due 
to the small fractional variability at these scales. Regularized 
linear inversion deals with this by smoothing the solution at 
the smallest time scales, resulting in a solution dominated by 
the longer multi-day variability features. For all the objects 
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Epoch (days) Epoch (days) Time delay (days) 


Figure 1. The effect of using continuum models with varying lev¬ 
els of small scale structure. The upper panels show the CAR(l), 
a = 1.0 in Equation continuum model used for the results in 
this paper. The panels below show results for continuum models 
with increasing a , resulting in decreasing structure at short time 
scales, corresponding to progressively steeper, or more negative, 
power spectral density slopes. The left panels show the contin¬ 
uum light curve of Arp 151 from LAMP 2008 as grey error bars. 
In each case, the best fit (median and la percentiles) continuum 
model is shown as a blue band and a single realization of the best 
fit continuum model is shown as a solid black line. The middle 
panels show the Arp 151 H (3 emission line light curve (black er¬ 
ror bars) along with the fit from RLI (green line and band) using 
the corresponding continuum models on the left. The right pan¬ 
els show the resulting response functions from RLI. The excellent 
agreement between response functions derived using different con¬ 
tinuum models leads us to conclude that the level of short time 
scale structure in the continuum model does not affect the RLI 
results derived from the LAMP 2008 dataset. 


analysed, the time scale of interest (a few days) is well within 
the sampling rate of the data. 

For each continuum light curve we find the best fit 
Gaussian proc ess parameters using the dnest3 Nested Sam¬ 
pling code bv iBrewer et al.l ( 2010 ). The best-fit Gaussian 
process is used to interpolate the continuum light curve, 
which we need to calculate the response for arbitrary time 
delays. 

In addition to interpolation, we use the statistical vari¬ 
ability of the Gaussian process to generate a number of re¬ 
alizations of the continuum light curve that are used to esti¬ 
mate statistical errors on the calculated response functions. 
This allows us to include measurement errors in the con¬ 
tinuum light curves that otherwise do not explicitly enter 
in the RLI formalism. For each continuum light curve, we 
generate 1000 realizations of the best-fit Gaussian process. 
This library of continuum realizations is used throughout the 
analysis to calculate response functions and response maps. 


3.3 Emission-line light curve errors 

Emission-line light curve errors cp(t) are included explicitly 
in the RLI formalism in eqs. © and ©• We allow for the 


possibility of additional systematic errors in the emission 
line fluxes by bootstrap resampling of the emission line light 
curves in each inversion. Bootstrap resampling is done by re¬ 
sampling data points in the light curve, such that each point 
can be chosen zero or more times and the total number of 
data points is kept constant. If a point is chosen N times, its 
error is reduced accordingly by y/N. We find that the results 
are only weakly affected by bootstrap resampling, indicating 
that systematic uncertainties affecting individual epochs in 
our data sample do not have strong influence on the results. 


3.4 Testing on simulated data (ID) 

We test our RLI method on simulated velocity-integrated 
(ID) light curves convolved with a selection of different 
response functions to produce a continuum-emission line 
light curve pair. To make the light curves as realistic as 
possible, the continuum light curve is simulated using a 
Gaussian process with a power spectral density similar to 
that found by recent Kepler observations (jEdelson et al.l 
120141 ). The normalization of the simulated continuum light 
curve is chosen to match t hat of typical LAMP 2008 light 
curves ([Bentz et al.l 1 2009c ), such that the simulated con- 
tinuu m has a fractional va riability of F var ~ 12 per cent 
(Fv ar ; IVaughan et al.l [20031 ) and a signal-to-noise ratio of 
SNR ~ 100. Emission line light curves are generated by 
convolving the simulated continuum light curve with the 
chosen response function (see Fig. El)- The simulated light 
curves are degraded and down-sampled to a level similar to 
the LAMP 2008 dataset (see section [2]) . We test our code 
on five different response functions: top-hat, multimodal 
top-hat, Gamma distribution, sinusoidal, and a delta func¬ 
tion. The simulated light curves and response functions are 
shown together with the fits from RLI in Fig. [2j We further 
test our method by simulating continuum light curves us¬ 
ing a damped random walk. Regularized inversion performs 
slightly better in this case, especially for small values of ft, 
due to the increased structure at short time scales, but the 
overall results are very similar to the simulated Kepler light 
curves shown in Fig. [2] 

RLI is generally successful in recovering ah simulated 
response functions, but like all reverberation mapping meth¬ 
ods the overall performance depends on the number of strong 
variability features in the light curve sample. Because of the 
smoothing constraint and the finite sampling of the light 
curves, RLI cannot fit very sharp features in the simulated 
response functions. This is particularly evident in the case of 
the top-hat and delta function response functions. For the 
smooth Gamma and sinusoidal distributions, we find that 
RLI is able to match the full shape of the response function. 
We tested the method on different simulations and found 
that the deviations at particular time delays are driven by 
noise in the input data and thus cannot be reduced without 
re-sampling, or extending, the light curve. 

The smallest scale resolvable by RLI can be seen in the 
recovery of the delta-function in the bottom panel in Fig. 
[2j This smallest scale, or point spread distribution, comes 
from the incomplete sampling of the light curves as well as 
the smoothing imposed on the solution to avoid fitting noise 
in the input. 

We test the effect of changing the regularization scale 
ft by one order of magnitude in each direction. The effect 
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of increasing the regularization (k/k o = 10, dashed blue 
line in Fig. is to smooth the response functions, thereby 
reducing some of the fluctuations in the wings, but also sig¬ 
nificantly increasing the width of the PSF (Fig. [2 lower 
right panel). In the same way, a reduction in regularization 
scale (k,/k>o = 0.1, dotted red line in Fig. [2} tends to amplify 
spurious features in the response functions while improving 
the PSF by decreasing its width. The best trade-off between 
high resolution and signal-driven results is a matter of pref¬ 
erence. We find that k — kq provides a good balance between 
resolution and noise. For this reason, and for consistency, we 
fix the regularization scale to n = ko for all our analysis. 

It is important to note that these simulations assume 
an exact linear relation between the continuum and emis¬ 
sion line, Equation ©. including only random errors and 
sampling gaps. The simulations and the results from RLI 
presented in Fig.[2thus represent an artificial best-case sce¬ 
nario. The main purpose of Fig. [2 is to illustrate the effect 
of changing the regularization scale k, show the finite point 
spread distribution of the RLI method due to discrete sam¬ 
pling, and the ability of RLI to recover negative and multi¬ 
modal response functions. 

3.5 Testing on simulated data (2D) 

In addition to testing one-dimensional response functions, 
we test our RLI implementation on data simulated us¬ 
ing the geometric and dynamical model of the BLR used 
i n dynamical modell i ng of reverberatio n mapping data 
([ Pancoast et al.l 1 20 111 : I Brewer et al.l 1201 lh . The dynamical 
modelling method simulates the broad line region as a num¬ 
ber of discrete point particles in a parametrized geometrical 
configuration. We test RLI on two different simulated BLR 
configurations, one in which the BLR particles are in near¬ 
circular orbits, and o ne with entirely inflowing orbits (see 
IPancoast et al ll2014al for details on the simulated datasets). 
The results are shown next to the true response maps in Fig. 
[3] We note that the data are simulated using a kinematic 
and geometric model for the BLR, and are not necessarily 
representative of true AGN response maps. 

These simulations provide a qualitative indication of the 
2D point-spread-function of the RLI method. Because the 
solution is forced to be smooth in the presence of noise, the 
resulting response maps will not reproduce the input maps 
exactly. Instead, we retrieve a smoothed version of the input 
response maps. The stripes in the RLI maps appear because 
we analyse each velocity bin individually. This means that 
the solution is not smoothed in the velocity (or wavelength) 
direction. We did this to simplify the calculations and to let 
any correlations between velocity bins be driven by the data 
only. This is contrary to maximum entropy, where the result 
is also smoothed (or regularized) in the velocity direction. 
If some velocity bins have exceptionally large errors, RLI 
might fail in that bin and produce a flat response. 

Although the response maps from dynamical modelling 
and RLI look somewhat different, they share some of the 
same symmetries. The effect of in-fall in the upper panel 
in Fig. [3] is clearly recovered in the RLI map. Likewise, the 
lower panel in Fig. [2 shows that dynamical modelling and 
RLI both show a symmetric response on either side of the 
emission-line centre, suggesting a structure where the BLR 
clouds are on closed orbits around the central black hole. 



Time (days) Time Delay (days) 


Figure 2. Testing regularized linear inversion on simulated data 
by recovering response functions for simulated light curves. The 
top left panel shows the simulated continuum light curve, gener¬ 
ated using a Gaussian process with a power spectr al density sim¬ 
ilar t o that found by recent Kepler observations (lEdelson et al.l 
l2Q14h . The magenta band plotted on top of the continuum shows 
the standard deviation of the Gaussian process fit used in the 
analysis. The left panels, below the continuum, show the emis¬ 
sion line light curves obtained by convolving the continuum with 
the corresponding response functions shown in the right panels. 
Before being analysed, the simulated light curves are degraded 
and down sampled to a level similar to the light curves of LAMP 
2008. The response functions recovered from regularized linear in¬ 
version are plotted on top of the input response functions in the 
right panels. The thin black line is the input response function 
while the solid green (red dotted and blue dashed) lines shows 
the response functions calculated from regularized linear inver¬ 
sion with a regularization scale of k/ko = 1 (k/ko = 0.1 and 
k/ko = 10). The response functions are, from top to bottom: top- 
hat, multimodal top-hat, Gamma distribution, sinusoidal, and the 
delta function. The grey dashed horizontal lines indicate the level 
of zero response. 

3.6 Effect of changing the regularization scale n 

Fig. H shows the effect of changing the regularization scale 
k, using the H ft data for Arp 151 as an example (see full 
results in Section l4.2.1D . We find that, for low vales of k, (i.e. 
more weight on the x 2 -term), the solution is very unstable 
to perturbations in the input, and the error in the response 
function increases at all time delays. The reason that the 
error is not a function of time delay, is that the noise is 
not correlated with the signal in the light curve. Therefore, 
the response function can draw power at any time delay 
to fit the noise, so long as it draws similarly less power at 
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Velocity (10 3 km s 3 ) Velocity (10 3 km s 3 ) 


Figure 3. Testing regularized linear inversion on simulated two- 
dimensional response map s (left column) generate d using the dy¬ 
namical modelling code of lPancoast et all (12014a! ). The response 
functions are calculated for each velocity bin individually using 
RLI, and then plotted together to produce velocity-delay maps 
(right column). 



other scales. Because the noise is random in each continuum 
realization, the time delays at which the response function 
draws its power will be evenly distributed in delay space, 
and the response function will fluctuate with more or less the 
same amplitude at all time delays. There will still on average 
be more power at scales correlated with the signal, thus the 
median shape of the response function peaks around the 
true delay, but the amplitude variations due to noise produce 
large error bars at all time delays. This behaviour is expected 
and illustrates the need for regularization to achieve stable 
solutions. 

As the regularization scale is increased, more weight is 
put on the smoothness of the solution. This suppresses the 
effects of noise in the input and emphasizes the large scale 
behaviour of the light curves. Because of the smoothing in¬ 
troduced in the regularization, there will be a minimum 
resolvable scale, analogous to an extended point-spread- 
function (see Fig. El. Any structure below the noise level 
will thus be smoothed out in RLI. 


4 RESULTS 

Here we present results from regularized linear inversion of 
light curves of five local AGN from the LAMP 2008 dataset. 
Fig. [5] to [IT] show results from regularized linear inversion 
of broad band photometric continuum light curves and H 
light curves from spectral decomposition. We provide fits 
to the light curves as well as integrated response functions 
and velocity-delay maps for all object. Table [I] lists derived 
median time delays for the integrated light curves of all AGN 


Figure 4. Effect of changing the regularization scale /■€, increasing 
the smoothness of the solution. Using Arp 151, H (3 as an example. 
The left panels show the H (3 light curve (black points with error 
bars) along with the fit from RLI (solid green line and band). 
The right panels show the corresponding response function found 
using RLI. The regularization scales affects all time scales in the 
response function equally, which is why the error bars are more 
of less constant as a function of time delay for each regularization 
scale k. 


analysed, as well as time delays from cross-correlation and 
javelin for comparison. 


4.1 Analysis 

4.. 1.1 Light curves 

The upper left panels in Fig. [5] to [lT] show the observed light 
curves plotted as black points with error-bars. The blue band 
on top of the continuum light curve shows the median and 
la percentiles of 1000 realizations of the Gaussian process 
best-fit model. This band indicates the range of continuum 
fight curve realizations used together with bootstrap resam¬ 
pling to determine errors on the derived response functions 
(see Sections E2 and 1531) . The middle left panels show the 
integrated emission-line light-curves obtained by integrating 
the corresponding emission-line over the wavelength range 
listed in Table [T] The solid line and green band on top of 
the emission line fight curves show the median and la per¬ 
centiles of the fits obtained from RLI when analysing all 1000 
realizations of the continuum together with 100 bootstrap 
resamplings for each continuum realization. 
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4.1.2 Integrated response functions 


4.1.5 Comparing with other methods 


The lower left panels of Fig. [5] to HU show the best-fit in¬ 
tegrated response functions under the smoothing constraint 
described in Section [3j The response function and the as¬ 
sociated error-bars are calculated from the median and la 
percentiles of RLI using 1000 iterations of continuum light 
curve interpolations as well as 100 bootstrap resamplings of 
the emission line light curve. Because we minimize not just 
the x 2 , but also the first derivative of the response function, 
the medians and error-bars of the response functions will be 
naturally correlated. 


4-1.3 Velocity-resolved response functions 

The upper right panels of Fig. [5] to on show the velocity- 
delay maps from regularized linear inversion of the spectra 
for all epochs in the data. Each spectral bin is treated indi¬ 
vidually, such that any correlations below the pixel resolu¬ 
tion (see Section lT2l) are driven by the data. As in the case 
of the integrated analysis, we fix the regularization scale to 
k, = Ko for all velocity bins. 

The velocity-delay maps show the response on a linear 
scale as a function of time delay and Doppler velocity with 
respect to the emission-line centre. White colour corresponds 
to zero response while black is maximum response (see the 
colour bars for each individual object). Below the response 
maps we plot the spectrum of the corresponding emission 
line. The black line is the mean spectrum across all epochs 
considered and the grey line is the root-mean-square vari¬ 
ability (standard deviation about the mean) of the spectra, 
normalized to the same scale as the mean spectrum. 

The lower right panels of Fig. l5l to fill show a selection 
of response functions for three separate velocity-bins. The 
location of the bins are indicated by the vertical dashed lines 
in the response maps and spectra. The colours correspond to 
the relative wavelengths, such that the red response function 
is for the highest radial velocity bin (redshifted with respect 
to the line centre), the green response function is for the line 
centre at 0 km/s and the blue response function is for the 
lowest radial velocity bin (blueshifted with respect to the 
line centre). 


4.1.4 The time delay 

The scalar time delay from cross correlation methods, tccf, 
is typically calculated as the centroid of the CCF above a 
threshold (usually 80 per cent of the maximum of the CCF). 
This effectively ignores any contributions from negative re¬ 
sponse. To compare our results to those from cross correla¬ 
tion, we calculate a RLI time delay, trli, as the median of 
the positive values of the response function. The time delays 
and error bars for each emission line and object, quoted in 
Table [T] and in the text, are then estimated as the median 
and la percentiles of all time delays determined from the re¬ 
sponse functions of all continuum realizations (see Sections 
I3.2l and l3.3l) . All time delays quoted in the text and in Table 
[l] are in the rest frame of the AGN. The time delay units on 
the axes of the figures are in the observed frame. 


We compare our results to those of other methods, including 
cross-correlation, javelin , maximum entropy, and dynam¬ 
ical modelling. 

The cross-correlation method calculates the cross¬ 
correlation function (CCF) between the continuum and 
emission line light curve to find the time delay, generally 
characterized by the centroid of the CC F, of the respond ¬ 
ing ^ as in the BLR. First applied by iGaskell fe Soarkel 


(j 19861 b t he method has si nce b een substantially devel¬ 
oped (e.g. Edelson &; Krolil^ 1 19881 : IWhite &; Petersonl Il994l ; 
IPeterson et al.1l998al b The CCF is related to the transfer 
function throu gh the light curve auto-correlation function 
(|Penstonlll99ll b This means that the width of the CCF de¬ 
pends on the variability characteristics of the AGN as well as 
the sampling of the data. For this reason, the CCF is rarely 
interpreted by itself, but is rather used as a tool to estimate 
the characteristic size of the BLR by calculating the time 
delay. Table [I] compares integrated time delays from RLI to 
time delays determ ined using c ross correlat i on met hods by 
iBentz etahl (|2009ch for H (3 and lBentz et al.l (|2010ah for H a 
and H 7 . Fig. [12] compares velocity-resolved time delays from 
RLI to velocity-r esolved time d elays from cross correlation 
of the H ^ line by IBentz et al.l (l2009ch . The comparisons for 
each object are described in Section 14.21 See the previous 
section for a description of how the RLI time delays are 
estimated from the response functions. 

javelin finds the ti me delay b y usin g a top hat model 
for the transfer function (|Zu et aljfeoill b Light curve vari¬ 
ability is modelled using a CAR(l) process, which is the 
same process used to model the continuum variability in 
our RLI method. Like the cross-correlation method, javelin 
provides only a scalar time d elay as a measure of the BLR 
structure. iGrier et al.l ( 2013b ) determined time delays using 
javelin on the LAMP 2008 dataset. We list time delays 
from IGrier et al.l (|2013bh in Table [T| alongside results from 
cross-correlation and RLI for comparison. All j avelin time 
delays tjavelin quoted in the text are also from IGrier et"ai1 
(l2013bl b 

_ T he dy namical modelling method ([Pancoast et al.l 

1 20 111 . l2014ah implements a full three-dimensional model 
of the BLR. The BLR emission is modelled as coming 
from a number of point particles that are drawn from a 
parametrized phase space distribution and linearly repro¬ 
cess radiation from a central source. All measured parame¬ 
ters have a direct physical interpretation, which means that 
the black hole mass Mbh and the virial factor / can be de¬ 
termined directly. Dynamical modelling results include fully 
velocity-resolved transfer functions. Because the method is 
based on a physical model of the BLR, the response will 
be correlated across velocity bins. This is contrary to our 
implementation of regularized linear inversion, where each 
velocity bin is analysed individually. 

Another important difference between dynamical mod¬ 
elling and RLI is that, because dynamical modelling traces 
photons through the BLR to the observer, the resulting 
velocity-delay maps are transfer functions representing the 
absolute reprocessed emission. This is contrary to RLI which 
subtracts the mean component of the light curves to only 
calculate the response in the emission-line to variations in 
the continuum, thus measuring response (rather than trans- 
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fer) functions. However, in the current imple mentation of 
dynamical modelling bv IPancoast et all (12011 ). the respon- 
sivity of the BLR gas is assumed to be constant throughout 
the BLR. In other words, the response of the BLR is corre¬ 
lated one-to-one with the emissivity distribution, regardless 
of the level of ionizing continuum. Because of this, the trans¬ 
fer functions and response functions will be directly propor¬ 
tional, and we can compare response maps from RLI to the 
velocity-resolved transfer functions of dynamical modelling 
directly. 

Maximum entropy methods (ISkill ing &; B rvanl Il984l : 
iKrolik et al.l Il99ll : iHorne et al.l Il99ll : iHornel Il994l ) work by 
selecting solutions to the convolution problem (Equation [2} 
that simultaneously provide a good fit to the data while be¬ 
ing as simple (smooth) as possible. The transfer functions 
are proposed based on a number of assumptions about the 
form of the solutions. This means that, like the dynamical 
modelling method, the derived response values will be de¬ 
pendent across velocity bins and the selection of the solution 
with the maximum entropy tends to result in very smooth 
velocity-delay maps. This is similar to RLI except that we 
impose no correlations between velocity bins, so any corre¬ 
lations between bins will be driven purely by the data. We 
compare our results to t he response map for H /3 in Arp 151 
from maxi mum entropy dBentz et al.ll2010bh and dynamical 
modelling ([Pancoast et alJ l2014bh in Fig. [T3j The compar¬ 
isons are described in Section [4.2.11 

4.2 Results for each AGN 

4-2.1 Arp 151 

Results for Arp 151 are presented in Fig. [5] for H/3 and Fig. 
[ 6 ] and 0 for H a and H 7 respectively. Light curves for Arp 
151 are the most variable and highest quality of the LAMP 
2008 dataset, which is why we include only Ho and H 7 for 
this object. RLI provides decent fits to the Arp 151 light 
curves, with some notable outliers at later epochs (e.g. H /3, 
Fig. [5]). Because of the smoothing constraint imposed on 
the solution we do not expect our code to fit these points. 
This is si milar to what is found in other analyses of the 
same data (|Bentz et al.ll 2010 bl : IPancoast et al.ll2014bh . Such 
strong variability on short time-scales is in any case not 
consistent with models where the BLR has only extended 
emission (see discussion in Section 15.ID . 

The integrated response function for H f3 (Fig. [5]) has 
a plateau from 0 days out to about 7 days where it 
starts to decrease. The shape of the response function is 
broa der than what is f ound by maximum entropy meth¬ 
ods dBentz et al.ll2010bh . We find a median H (3 time delay 
of trli = 4.0!q's days, consistent with results from cross¬ 
correlation (tccf = 4.01q 7 days) and javelin (tjavelin = 
3 - 6 io '.2 days). The integrated response function for H a (Fig. 
© has a flat low response out to about 6 days after which 
is rises to a peak response at a time delay of around 8 days, 
and then drops to slightly negative response after 10 days. 
The time delay for H a (trli = 6 . 8 ^ 4 ) is consistent with 
that from cross correlation (tccf = 7.81lq). The integrated 
response function for H 7 (Fig. [TJ) shows significant response 
at zero delay with a slight rise to a peak at 3 days after which 
the response drops to zero beyond 8 days. The time delay 


from RLI (trli = 3.01q'J) is somewhat larger, but consistent 
with, the result from cross correlation (tccf = 1.41q;®). 

The velocity-delay map for H f3 (Fig. [SJ) shows that the 
bulk of the response is centred on the emission line with 
a time delay of about 5 days, similar to the median de¬ 
lay calculated from the integrated response. There seems to 
be an area of increased response redward of the line centre 
from 0 — 1000 km/s. The velocity-delay map for H a (Fig. 
© shows a strong response around 10 days centred on the 
emission line and extending to lower time delays. The origin 
of the prompt response redward of the line centre is un¬ 
clear. It may be a spurious feature due to numerical effects, 
a poorly subtracted continuum (the continuum for H a was 
subtracted using a linear fit rather than spectral decompo¬ 
sition as for H /?, see Section lT2l) . or perhaps due to residual 
contamination from [N 11 ], A6586 A. We do not find evidence 
for a blueward plume from 15 to 20 days (not shown) a s 
seen in the maximum entropy maps in lBentz et al .1 (2010b). 
The velocity-delay map for H 7 (Fig. 0) shows a broad re¬ 
sponse with the longest delays at the centre of the line, and 
progressively decreasing time delays in the wings. We find 
additional prompt response at the centre of H 7 . 

The decrease of the median time delay in Arp 151 
from H a to H (3 a nd H 7 has bee n previously observed 
(|Bentz et alJl 2010 bl . see lGaskefill2009l for a review), and may 
be an effect of the varying optical depth for the Balmer lines. 
Because the optical depth is larger for the transitions be¬ 
tween lower excitation states, H a photons will have a harder 
time escaping the BLR without being absorbed. This results 
in the H a emission being predominantly directed back to¬ 
wards the ionizing source at the centre, whereas the lines 
formed by the higher excitation states, H /3 and H 7 , are pro¬ 
gressively more isotropic. That is one possible mechanism by 
which the median time delay can decline in higher excita¬ 
tion emission lines, while all the Balmer lines originate from 
hydrogen gas located at similar distances from the ionizing 
source. 

Fig. 1121 shows a comparison between RLI and velocity - 
resolved cross-correlation time delays bv lBentz et al.l (|2009cJ j 
for a number of velocity bins across the H /3 emission line. 
We find excellent agreement with cross-correlation for all 
velocity bins and confirm the signature of prompt response 
redward of the H /3 line centre, while the bulk of the response 
is at the centre of the emission fine at a time delay of around 
5 days. 

Fig. [I3l shows velocity-delay maps for H /3 from dynam¬ 
ical modelling, maximum entropy, and RLI. All methods 
find prompt response in the red wing of the H /3 emission 
line. This is a characteristic feature of BLR models with 
free-falling gas or a disk of gas contain ing a hot spot (see 
I Welsh fe Horneiri99ll : iBentz et alJ[ 2010 b ). Some models also 
produce prompt red-side resp onse for outflowing g as at par¬ 
ticular observed orientations (|Bottorff et alJ fl997h . The re¬ 
sult from RLI has a stronger response at the line centre, and 
a weaker prompt response in the red wing compared to dy¬ 
namical modelling and maximum entropy. In addition RLI 
finds prompt response at the line centre, which is not seen 
in the maximum entropy maps. 
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Velocity (10 3 km s x ) 




Figure 5. Results from regularized linear inversion of light curves from LAMP 2008 of the Arp 151 H f3 emission line. The left panel 
shows the continuum (upper left panel, black points) and integrated emission-line light curves (middle left panel, black points) along with 
the integrated response function (lower left panel, black points). The blue shaded region on top of the continuum light curve shows the 
median and lcr percentiles from the Gaussian process realizations used to model the uncertainty in the data. The green line on top of the 
emission line light curve shows the best fitting result from RLI, with the green shaded band indicating the lcr percentiles. The horizontal 
grey dashed line in the lower left panel indicates the location of zero response. The right panel shows the velocity-delay map (upper 
right panel) calculated by regularized linear inversion of light curves for each wavelength bin individually. Below the velocity-delay map, 
the mean (black line) and RMS (grey line) spectrum is shown as a function of velocity with respect to the H (3 emission line centre 
(middle-right panel). Normalized response functions for three velocity bins are shown below the velocity—delay map (lower right panel) 
with the bins indicated by dashed vertical lines in the velocity-delay map. The colours (blue, green, red) correspond to the Doppler shift 
with respect to the observer, with blue being negative velocity, green is the central velocity, and red is positive velocity. All time delays 
in the results figures are in the observed frame. 
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Figure 6. Same as Fig. [5]but for Arp 151, Ha. 
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Object 

Continuum 

Band 

Emission 

Line 

Wavelength 
Range (A) 

tccf 

Time delay 
(days) 
tjavelin 

trli 

Arp 151 

Arp 151 

Arp 151 

Mrk 1310 

NGC 5548 

NGC 6814 

SBS 1116+583A 

B 

B 

B 

B 

V 

V 

B 

Ha 

H/3 

h 7 

H/3 

H/3 

H/3 

H/3 

6575 - 6825 

4792 - 4934 

4310 - 4393 

4815 - 4914 

4706 - 5041 

4776 - 4936 

4797 - 4926 

7.8tfS 
4-o ±8;? 

1 4+ 0 - 8 

1 -0.7 

q y+0.6 

4-2±?;| 

2.3+j? 

3- 6+1 

4- 2+1 

5.5+;® 

7 4+ 0 - 1 

1 -0.1 

2 4+0.9 

Z -0.9 

6.8/+I 

4-0+s 

3-0+s 

O 7+0.3 

1 -0.3 

4 7+ 1 - 8 

1 -1.8 

7 Q + l- 2 

1 ■ c> —1.0 

2 fC 1 ' 1 
Z ' U -0.6 


Table 1 . Time delay s from cross-correl ation (tccf) are reproduced from iBentz et all (l2010alb Time delays from javelin (tjavelin) 
are reproduced from lGrier et al.l (2013b). All time delays are given in the rest frame. 




Velocity (10 3 km s *) 




Time Delay (days) Time Delay (days) 


Figure 7 . Same as Fig. [5]but for Arp 151, Hy. 


4.2.2 Mrk 1310 

Results for Mrk 1310 are presented in Fig. [8] There is more 
scatter in the light curves of Mrk 1310 compared to Arp 151. 
Because of the smoothing constraint, RLI is not able to fit 
the emission line light curve exactly. This is to be expected, 
as the method mainly picks up large scale variability to avoid 
problems with noise on shorter time-scales. The scatter in 
the H p light curve towards the end of the campaign is in any 
case difficult to reconcile with a simple linear response to the 
continuum light curve variations near the same epochs. 

The integrated response function for H (3 is more peaked 
than that of Arp 151. The scatter in the derived response 
functions from our error analysis is also smaller, indicat¬ 
ing that the inversion is stable under perturbations in 
the input light curves. The median time delay of the re¬ 
sponse function is trli = 2.71q'J days. This delay is slightly 
shorter than, but consistent with, the cross-correlation re¬ 
sult (tccf = 3.71q;6 days), while the result from javelin 
(tjavelin = 4.2+q i days) is longer by 1 day compared to 
the time delay from RLI. The integrated response function 
dips slightly below zero at longer time delays, which could 
indicate a negative response in the emission-line, although 


we do not find this to be significant (see discussion in Section 

EU). 

The velocity-resolved response map shows a strong re¬ 
sponse at the line centre with a time delay around 3 days, 
similar to the integrated time delay. The width of the re¬ 
sponse is about ±1000 kms -1 , similar to Arp 151. The delay 
is fairly constant across the H ft line. This result agrees well 
with velocity-resolved cross-correlation that shows a nearly 
flat response as a function of velocity, with slightly shorter 
delays in the wings at ±1000kms -1 (Fig. fl2l) . It is also 
consistent with dynamical modelling results, with a peak 
response centred on the emission line at fairly short time 
delays (Fig. [13l middle left panel). 

4.2.3 NGC 5548 

Results for NGC 5548 are presented in Fig. [9] The H light 
curve of NGC 5548 has more scatter relative to variabil¬ 
ity amplitude compared to Arp 151 and Mrk 1310. Conse¬ 
quently, RLI has a more difficult time fitting this object. 

The integrated response function increases from zero 
delay up to around 3 days and then has a slowly decreasing 
plateau out to about 8 days after which it drops off. This 
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Time Delay (days) 


Time Delay (days) 


Figure 8. Same as Fig. [5]but for Mrk 1310, H/3. 


yields a median time delay of trli as* 4.7l J'J days, consistent 
within la with cross-correlation (tccf = 4.2+^g days) and 
JAVELIN (tjavelin 5.5to® days). 

The H/3 emission line in NGC 5548 is significantly 
broader compared to the other objects analysed and the 
variability across the line shows more irregular structure, as 
is seen in the velocity-resolved response map in the upper 
right panel of Fig. [9] We find a somewhat weak response 
at the line centre, which peaks at around 5 days. On ei¬ 
ther side of the line centre at ±5000kms _1 we find stronger 
isolated response features with peak time delays around 8 
days. The velocity-resolved time delays agree well with cross¬ 
correlation (Fig. fl2l) . Dynamical modelling and RLI both 
show a level of prompt response in the red wing of H/3, 
while the rest of the maps show little similarity fFig. H3l) . 


4.2.4 NGC 6814 

Results for NGC 6814 are presented in Fig. [TOj The light 
curves of NGC 6814 show an appreciable amount of vari¬ 
ability. The continuum R-band light curve has a clear broad 
peak at HJD — 2454000 = 560 days after which it drops 
off to a slowly rising plateau that extends to shortly after 
HJD — 2454000 = 600 days. The emission line light curve 
has the same over-all trend as the continuum, but instead 
of having a plateau after the initial peak, it rises to almost 
the same level as the initial peak. Because of this, RLI has 
a difficult time matching the first and second peaks in the 
emission line light curve, explaining why the fit is underes¬ 
timated for the second peak. 

The integrated response function peaks around 8 days 
and has broad flat wings to either side (Fig. 1101 lower 
left panel). The median time delay of trli = 7.31J q days 
is consistent with results from cross-correlation (tccf = 
6-5l?'.o days) and javelin (tjavelin = 7.4^;} days). 

The velocity-resolved response map for NGC 6814 
shows a clear isolated response peak centred on the H/3 


emission line (Fig. [TUI upper right panel). There is a some¬ 
what peculiar dip in the response map just blueward (~ 
—100 km s -1 ) of the line centre. Further away at around 
±500 km s -1 the response peaks on either side of the line 
centre and then drops off to zero beyond ±1500kms _1 . 
The velocity-resolved time delays agree well with cross¬ 
correlation at centre of the H/3 line (Fig. 1121 ) . At lower ve¬ 
locities we find a slightly longer time delay compared to the 
cross-correlation result. At higher velocities, blueward of the 
line, RLI was unable to recover a time delay due to a lack of 
emission line response. Looking at the fully resolved response 
map from RLI in Fig. [TU]we see that RLI calculates little 
response in the wings beyond ±1500 km s -1 , which might 
explain the discrepancy with the cross-correlation results. 
Dynamical modelling of NGC 6814 tends to prefer shorter 
time delays than what is found using other methods (see 
Fig. [T3l) . The highest velocity bin failed to produce a time 
delay estimate due to a noisy response function. 


4.2.5 SBS 1116+583A 

Results for SBS 1116+583A are presented in Fig. [TT] The 
emission line light curve of SBS 1116+583A lacks significant 
features above the noise level. As a result the RLI fit does 
not match the emission line light curve well. The results for 
this object, including the velocity-delay map, should thus 
be interpreted with caution. 

The integrated response function prefers a relatively 
short time delay of trli = 2.0^ J'g days, which is slightly 
shorter, but fully consistent with cross-correlation (tccf = 
2.3^0 5 days) and javelin (tjavelin = 2.41 q '.9 days). 

The velocity-resolved response map for SBS 1116+583A 
(Fig. 1111) shows a multimodal response map with a broad 
response component at short time delays, around 2 days, 
and an additional isolated response at 10 days at the centre 
of the H /3 line. This is a good example of the flexibility of 
RLI to reconstruct complicated velocity-delay maps. The 
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Figure 9. Same as Fig. [5]but for NGC 5548, H (3. 
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Figure 10. Same as Fig. [5]but for NGC 6814, H/3. 


velocity-delay map shows some resemblance to models in 
which the BLR is confined to an iso tropically illuminated 
disk (see Fig. 5 in lBentz et al.1 201 Obi ). 

It is not very meaningful to calculate a single time delay 
from a multimodal response function, even so we apply the 
same procedure for determining the time delay as for the 
previous objects (see Section |4. 1.41) . The resulting velocity- 
resolved time delays agree well with cross-correlation, show¬ 
ing longer time delays at the centre of the emission line and 
shorter time delays in the wings fFig. fl2l). We note that the 
time delay calculated close to the H (3 line centre fall between 
the two modes in the response map, and is thus not well con¬ 
strained. We include these time delays for comparison with 


the velocity-resolved cross-correlation, but they should be 
interpreted with caution, as suggested by the large error- 
bars. The lowest velocity bin failed to produce a time delay 
estimate due to a noisy response function. 


The results for SBS 1116+583A from RLI are in general 
agreement with dynamical modelling, although dynamical 
modelling does not reproduce the same strong response at 
t = 10 days and generally prefers shorter time delays at the 
line centre (see Fig. fl3l). 
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Time Delay (days) 


Figure 11. Same as Fig. [5]but for SBS 1116+583A, H (3. 
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Figure 12. Comparison of H /3 time delays from regularized linear inversion with cross-correlation results from lBentz et all ll2009ch . The 
upper panels show time delays as a function of velocity with respect to the H (3 emission line ce ntre. Red filled circles with error bars show 
time delays from this paper. Black open squares with error bars are from lBentz et all (|2009cT) . Light curves are obtained by integrating 
the emission in each velocity bin indicated by the horizontal error-bars. The time delays from RLI are calculated as the median of the 
positive part of the response (see Section [4.1 Ah . The lower panels show the mean (black) and RMS (grey) spectra. All time delays in 
this figure are in the observed frame. 


5 DISCUSSION 


5.1 Assumption of a constant linear response 

Long temporal baseline multi-epoch observations of AGN 
reveal significan t changes in the broad emission line equiva¬ 
lent widths (e.g. Kinney et al.lll9 90l: Pogge Petersonlfl992l : 


iGilbert &; Petersonll2003l: iGoad et al.ll2004h and line profiles 

(e.g. IWanders &; Petersonl Il996l : ISergeev et al.1 l200ll . 120071 ) 
on time-scales of months to years. These effects could possi¬ 
bly account for some of the discrepancies found when match¬ 
ing longer time-scale variability patterns using linear rever¬ 
beration mapping, such as the inability of RLI to match the 
second peak in the emission line light curve in NGC 6814 

(Fig. [B. 


In addition to long time-scale non-linearities, some light 
curves of the LAMP 2008 dataset show high-variability fea¬ 
tures (e.g. late epochs in Mrk 1310, Fig. if these fea¬ 
tures are associated with the BLR, they are inconsistent 
with models in which the BLR emission comes from an ex¬ 
tended region. This remains true even after exclusion of po¬ 
tentially unreliable spectra in the LAMP 2008 dataset, as 
described in section EZ2 Whether these outliers are due to 
systematic measurement uncertainties affecting individual 
epochs, an indication of unknown processes in the BLR, or 
a combination of these, is unclear. 

Whatever the origin of non-linear features in the emis¬ 
sion line light curves, current reverberation mapping tech¬ 
niques, such as cross-correlation, maximum entropy, and 
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Figure 13. Comparison of velocity-dela y maps from regularize d linear inversion to maximum likelihood for Arp 151 (iBentz et al.ll2010tJ ) 
and dynamical modelling for all objects (jPancoast et al.ll2014bl l. All three methods used the same dataset from LAMP 2008 (see Section 
0. Regularized linear inversion and dynamical modelling used spectral decomposition for continuum subtraction, while the maximum 
likelihood analysis used a linear fit to remove the continuum level below the line. 


RLI, will be insufficient to describe them. Even so these 
methods remain valuable tools for testing models of the 
BLR and investigate departures from linearity in AGN light 
curves. Implementation of photoionization physics into dy¬ 
namical modelling codes may be a way to probe non-linear 
processes in the BLR. 

5.2 Ionizing versus observed continuum 

Studies suggest a non-negligible time delay between the ion¬ 
izing UV continuum and the optical cont inuum on the order 
of 1 day in NGC 7469 and NGC 5548 (e.g.lCollier et al.l[l998l : 
iMcHardv et al.ll20l4lEdelson et al.ll2015h . This means that 
the V and B band continuum light curves originate from a 
region comparable in size to some of the shorter time de¬ 
lays found for Balmer emission lines. This extended contin¬ 
uum emission region will introduce geometric smoothing in 
the optical continuum light curves used for reverberation 
mapping. If the delay between the ionizing continuum and 
the measured optical continuum is constant on time-scales 
of reverberation mapping campaigns, linear reverberation 
methods can still be applied, but the interpretation of the 
response functions will have to be revised to include not 
just the geometry of the BLR, but also the geometry of the 


extended optical continuum emission region. It will be in¬ 
teresting to extend this analysis in the future by performing 
RLI of optical emission line variability as driven by UV con¬ 
tinuum to assess the magnitude of these systematic effects. 

5.3 Response function errors 

The smoothing condition used to regularize the inversion 
problem introduces correlations in the response functions 
by linking nearby points through the first-order derivative. 
This correlation is an implicit assumption of the method 
and it reflects our belief that the emission line light curve is 
produced by a superposition of photons emitted from a fairly 
homogeneous distribution of gas in the BLR. The error bars 
on the response functions should thus not be interpreted as 
errors on the individual points of the response function, but 
rather as an indication of the range of solutions we get by 
re-interpolating the continuum light curve and re-sampling 
the emission line light curve. 

In our current implementation of RLI no correlation in 
response is imposed between response velocity bins. This 
can be witnessed by the vertical streaks in the presented 
response maps. Not imposing a correlation has the advan¬ 
tage that any correlation between velocity-bins in the final 






























Constraints on the BLR from regularized inversion 17 


response maps must be driven by the data, or at least sys¬ 
tematic effects correlated with the emission lines. It seems 
natural to assume that the BLR emissivity is correlated in 
position as well as velocity space, motivating a smoothing 
constraint in the velocity domain, as is also assumed in 
dynamical modelling and maximum entropy methods. Be¬ 
cause we found the velocity-bins to be naturally correlated 
in the response maps, and in order to keep the method fast 
and simple, we decided not to impose any smoothing in the 
velocity-domain. 

5.4 Negative response values 

Some models suggest that negative values occur naturally 
in AGN response func tions (|Sparkelll993l : IGoad et al.lIl993l : 
iGoad &; Koristal l2014h . and negative response value s have 
recently been obs erved in X-rays (iMcHardv et al.l 120071 : 
iFabian et alJl2QQ9h . Negative response can arise in different 
ways, such as if an increase in ionizing flux causes part of the 
responding gas to become fully ionized, or if the BLR struc¬ 
ture changes on a dynamical time-scale, such as if the contin¬ 
uum decreases while new material falls into the BLR. Note 
that in both of these examples the BLR structure changes on 
time-scales of the same order as the time delay. This means 
that the assumption of a constant linear response breaks 
down and more general methods, such as dynamical mod¬ 
elling, are required to account for the time varying part of 
the response. 

An advantage of RLI over other reverberation mapping 
methods, is that it does not require the response function to 
be everywhere positive. We find some preference for nega¬ 
tive response in some of the objects analysed (e.g. Mrk 1310, 
Fig. [ 8 ]). While none of the response functions show strong 
evidence for negative response these results may be an indi¬ 
cation that more complicated processes are taking place in 
the BLR, than what can be contained in a simple constant 
linear response model. However, negative values occur natu¬ 
rally as artefacts when linear inversion methods are applied 
to discrete and noisy data, even if the true response is ev¬ 
erywhere positive. This is somewhat in analogy to aliasing 
distortions when dealing with discrete Fourier transforms. 
These artefacts are seen as ringing effects, such as observed 
in the inversion of the simulated ID data (see Fig. [2]). It is 
not possible to impose positivity while keeping an anal ytical 
expression for the response function (|Vio et al .11 19941 ). and 
since simplicity and flexibility were the main motivations for 
using RLI, we leave the positivity constraint to other meth¬ 
ods such as maximum entropy and dynamical modelling. 
With these considerations in mind we emphasize that all 
the negative response values found in our analysis of the 
LAMP 2008 data are of low statistical significance. In par¬ 
ticular, the integrals of the response functions are positive 
for all objects analysed. 


6 SUMMARY AND CONCLUSION 

Building on previous work bv lKrolik fe Donel (|l995h . we de¬ 
velop a new method for reverberation mapping based on reg¬ 
ularized linear inversion that includes statistical modelling 
of the AGN continuum light curves. In this implementation, 
response function errors are evaluated using a combination 


of Gaussian process fitting of the continuum light curves 
as well as bootstrap resampling of the emission line light 
curves. Regularized linear inversion has the advantage over 
other reverberation mapping methods that it makes few as¬ 
sumptions about the shape of the response function, and it 
allows for negative values in the response function. In addi¬ 
tion, because the method is based on an analytical solution 
to the transfer equation, it is very computationally efficient. 
The scale of regularization is a function of the data only, such 
that no user-input is required, except for the time baseline 
of the analysis. This means that flexible response functions 
can be calculated quickly and consistently. 

We test the method on simulated data and show that it 
is able to recover unimodal and multimodal response func¬ 
tions as well as response functions with negative values (see 
Fig. [2}. We further test the method on data simulated us¬ 
ing dynamical modelling and show that it is able to recover 
2-dimensional velocity-delay maps (see Fig. 0. 

We apply RLI to light curves of five nearby Seyfert 
1 galaxies, taken by the LAMP 2008 collaboration, and 
present time delays, integrated response functions and 
velocity-delay maps for the H f3 line in all objects, as well 
as Ha and Hy for Arp 151. This is the first time a reverber¬ 
ation mapping method allowing for negative response has 
been applied to a large dataset. Our results are in good 
agreement with previous studies based on cross-correlation 
and dynamical modelling, offering a powerful corroboration 
of the assumptions underlying these reverberation mapping 
methods. 

The main physical results of this work can be summa¬ 
rized as follows: 


(i) Despite using a very different method for calculating 
time delays we find that our results are in good agreement 
with results from cross correlation. 

(ii) We find asymmetric response of the H [3 emission line 
in Arp 151, with prompt response in the red wing of the 
emission line, consistent with models that include bulk gas 
flows in the BLR. 

(iii) We confirm previous studies finding that lines origi¬ 
nating from high excitation states, such as Hy, have shorter 
time delays compared to lines originating from lower excita¬ 
tion states, such as Ha and H /3. 

(iv) While some objects, such as Arp 151 and Mrk 1310, 
show a degree of negative response at longer time delays, we 
find no conclusive evidence for negative response values in 
any of the objects. 


This work demonstrates that regularized linear in¬ 
version is a valuable tool for reverberation mapping and 
a worthwhile complementary method for analysing high 
quality reverberation mapping datasets. The recent multi¬ 
wavelength reverberation mapping campa ign of NGC 5548 
(|De Rosa et al . 20151 : lEdelson et al.| [i2015j. including simul¬ 
taneous ground-based monitoring in the optical and spec¬ 
troscopy in the ultraviolet with the Cosmic Origins Spec¬ 
trograph on the Hubble Space Telescope , will provide excep¬ 
tionally high quality data and thereby a great opportunity 
to test the RLI method and its underlying assumptions. 
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